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Abstract 

We present an updated version of the SimProp Monte Carlo code: a 
simulation scheme to study the propagation of ultra-high-energy cosmic 
rays through diffuse extragalactic background radiation. The new version 
of the code presents two important updates: (i) it treats in a full stochastic 
approach all interaction channels involving ultra-high-energy cosmic rays 
and (ii) it takes into account the production of secondary (cosmogenic) 
neutrinos. This new version of SimProp was tested against different sim¬ 
ulations code, in particular the production of secondary neutrinos was 
compared with the fluxes expected in the scenario of Engel, Seckel and 
Stanev. 
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1 Introduction 

The present technical note should be intended as a user-guide that discusses 
modifications and improvements of the SimProp MC code respect to its first 
released version, SimProp v2r0, introduced in Ref. pQ. The code implemented 
here, SimProp v2r2, presents several changes regarding the propagation of both 
UHE extragalactic protons and nuclei, with a complete new module devoted to 
the computation of secondary neutrinos produced during the propagation. An 
intermediate version, SimProp v2rl, was used in Ref. |2] and briefly described 
there. Compared to it, SimProp v2r2 fixes few bugs affecting execution times 
and the low-energy tail of neutrino spectra and implements, in a full stochastic 
approach, the photo-pion production process (for both protons and nuclei) on 
the extragalactic background light (EBL), previously neglected, and on the cos¬ 
mic microwave background (CMB). SimProp is a free open-source code, available 
upon request by writing to SimProp-dev@aquila.infn.it 

A SimProp run consists of N events (for a user-supplied value of N ), each 
consisting in the generation of a primary proton or nucleus and its propagation 
from the source to the observer along with that of any secondary particles 
produced. All particles are assumed to travel rectilinearly (without taking into 
account the possible presence of intergalactic magnetic fields) at the speed of 
light, so only one coordinate, the redshift z, is used to keep track of time and 
positions. Each entry of the output file lists the initial energy, redshift, mass 
number and atomic number of the primary and an array with the final energies 
and types of the particles reaching the observer. 
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2 The command-line input 


SimProp v2r2 recognizes the following command-line parameters: 

-s Seed of the random number generator (default: 65539). 

-N Number of primary protons or nuclei to be injected (default: 100). 

-L Choice of the EBL model. 0: none; 1: Stecker et al. 0] (default); 2: power- 
law approximation of the previous model, as in 00; 3: Kneiske et al. 

0. 

-A Mass number of primaries, A; n j (choosen at random with - A 0; default: 56). 

-S -1: pion production approximated as continuous for protons, neglected for 
nuclei (as in SimProp v2r0); 0: pion production approximated as continu¬ 
ous for both protons and nuclei; 1: pion production treated stochastically, 
on the CMB only (default, as in SimProp v2rl); 2: pion production treated 
stochastically on both CMB and EBL. 


-D 0: all nuclei with mass number A are treated as the corresponding beta- 
decay stable isobar (as in SimProp v2r0); 1: neutrons and unstable nuclei 
can also be produced, but they are assumed to immediately undergo beta 
decay (default). 

-e Minimum log 10 (-E/eV) at injection (default: 17). 

-E Maximum log 10 (E/eV) at injection (default: 21). 

-z Minimum injection redshift (default: 0). 

-Z Maximum injection redshift (default: 1). 

-o 0: output as in SimProp v2r0 (default); 1: compact output, described below 
in Sec. 0 2: both output formats. 


The initial redshift and logE of primary particles is always uniformly dis¬ 
tributed between the limits given via the -z, -Z, -e and -E switches. In order 
to study different distributions, each event has to be weighed by an appropriate 
function of £j n j and Z[ n j; for instance, weighing each event by the function 


■j(Ei 


E 


1-7 


mj 


mj) ~inj ) 


(1 + z)\/ (1 + z) 3 fl m + 12 a 


(1) 


corresponds to assuming identical sources, uniformly distributed in comoving 
volume, each emitting particles with a power-law spectrum dN oc E^dE- m y 


3 The algorithm 

During each event, a stack contains all the particles to be propagated. At 
the beginning of the event, the stack only contains the primary particle, with 
user-specified mass number and with initial energy and redshift sampled from 
the user-specified ranges. New particles are added to the stack when produced 
during the propagation, and particles that interact, decay, or reach Earth are 
removed from it. The event is over when the stack becomes empty. 
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3.1 Propagation of protons and stable nuclei 

When propagating a proton or a stable nucleus the redshift interval between its 
production point z pro d and 0 (Earth) is divided into steps zo = z pr od, zi, z n = 
0, shorter near the production point than near Earth. During each step, there 
are two types of processes the particle can undergo: 

• those which can be approximated using the continuous (deterministic) 
energy loss approximation EEI and do not involve the production of any 
new particle to be tracked, namely the adiabatic loss due to the expansion 
of the Universe (redshift loss) and the process of electron-positron pair 
production on the CMB photons; 

• those which are treated as discrete (stochastic) interactions, with the in¬ 
teraction point, energies and types and/or number of outgoing particles to 
be sampled stochastically, namely the processes of photo-pion production, 
involving protons and nuclei, and photo-disintegration, involving only nuc¬ 
lei. 


3.1.1 Continuous energy losses 

At each redshift step (zj_i, z t ], the continuous energy losses are simply treated 
by numerically integrating from Zj_i to Zi the differential equation for lnr 


din r 

dz 




( 2 ) 


where T is the Lorentz factor of the particle, ft is the fractional energy loss per 
unit time, and 

_ 1 (3) 

dz Hq( 1 + z)y/ (1 + z 3 )f! m + Ua 

fixes the cosmology. In our computation we use D m = 0.3, 12 a = 0.7, and Hq = 
7.11 x 10 _11 /year ss 70 km/s/Mpc.The function ft in ([2]) is the sum of two terms, 
one for the redshift loss and one for electron-positron pair photoproduction. The 
former is computed as 


/3 ad (z) = H(z) = H ov / ( l-(- z) 3 n m + U A , (4) 


and the latter is 


Z 2 

ft pair(Z, A, r, z) = — (1 + z) 3 /3 pair (proton, (1 + z)r, z = 0), (5) 

where /3 pa ir for protons at z = 0 is interpolated from a list of tabulated values, 
computed as described in Ref. 7]; for non-proton nuclei, Z 2 /A is approximated 
as A /4 (exact when A = 2 Z). 

3.1.2 Discrete interactions 

The following scheme is used to decide whether and when the particle undergoes 
a discrete interaction and, if it does, the type and the products of the interaction. 
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Sampling the interaction point At the beginning of the propagation, a 
random number u is sampled from the uniform distribution between 0 and 1. 
The probability that the particle survives at redshift z without interactions is 
given by 



( 6 ) 


where the total interaction rate (probability per unit time) 1/r is computed 
at each step z- t as described below, and the integral is approximated via the 
trapezoidal rule, i.e., 


-In ^ 


bipi_i 


1 

2 




{Zi - Zi- 1 ) . 


(7) 


If at the end of a step pi < u. the particle is considered to have interacted 
during that step; the interaction point Zi nt is found by linearly interpolating p 
between z,_i and Zi and solving for p(z- ln t) = u, and the interaction energy E- m t 
is found by integrating m from Zi-\ to Zi n t. These are used to sample the 
number, type and energy of the outgoing particles as described below, adding 
these particles to the stack. 

If at the end of the last step p n > u , the particle is considered to have 
reached Earth; its mass number, atomic number, and final energy are recorded 
in the output file. 


Interaction rate The total interaction rate is given by 
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1 
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de de', 


( 8 ) 


where cr(e') is the total cross section for interactions with photons with energy e! 
in the nucleus rest frame (NRF), e( nin is the lowest value of e' at which the 
interaction is kinematically possible, and n 7 (e) de is the number per unit volume 
of background photons with energy between e and e + de in the laboratory 
frame. (The photon energy in the NRF is related to that in the lab frame by 
e' = re(l — cos 6), where 6 is the angle between the nucleus and the photon in 
the NRF; therefore, 0 < e' < 2Te.) 

We compute © as the sum of a term for pion production on the CMB 
T pion cmB’ one f° r pi° n production on the EBL T^ on EBL , and one for photodisin¬ 
tegration tjA. We assume that a nucleus behaves as A independent nucleons in 
pion production, i.e., Tpio n (A T, z ) = j 4' r I )ion(P r0 ^ 0n ’ z), because the energies 
involved are much larger than the binding energy of nucleons. This assumption 
corresponds to neglecting the nucleus recoil BUS]. 

We introduce the quantities 


7(e) 



n 7 (e) 

2e 2 


de 


(9) 


and 



m 2 )cr(s') ds' = 4m 2 J e'a{e') de ', 


( 10 ) 






GSSI/PHYS/2015.0042 


6 


where m is the mass of the particle (a nucleus in the case of disintegration and 
a nucleon in the case of pion production) and s is the centre-of-mass energy 
squared s = m 2 + 2me'; m can be also written as 

k = /2 A {m2+4mrc) ^ de < n > 

= f/, (12) 

c min x ' 

The photon background ro 7 is the sum of two terms, one for the CMB and 
one for the EBL. The CMB spectrum is known exactly at all redshifts, being 
a black-body spectrum with temperature T = (1 + z)T$, where To = 2.725 K, 
corresponding to a photon density 

" cmb(e) = iexp( e /4r) - 1 <13) 


and 

kviT 

IcMB{e) = ln (! - exp(-e/fc B T)); (14) 

this implies that T^ niCMB (T, z) = (1 + z) 3 Tj2* niCMB ((l + z)T,z = 0). On the 
other hand, EBL is not precisely known, and needs to be approximated using 
phenomenological models; for this purpose, depending on the command-line 
settings, 1) /ebl(c, z) is interpolated from a 2D grid as a function of loge and 
z, obtained from numerically integrating ([5]) using values based on the model of 
Ref. [3], or 2) the same power-law approximation is used as in Ref. [21 E], or 3) 
^ebl(c,~ = 0) is interpolated as a function of e from the model of Ref. 6 ), and 
the interaction rates at z > 0 are approximated as those at z = 0 times a scale 
factor interpolated as a function of z. 

We used the photo-pion production cross section for protons er p i on (e') as 
computed by the SOPHIA Monte Carlo code [5j. Performing a numerical integ¬ 
ration of m we built a table of values from which < f , p i on (s) can be interpolated. 
The photo-disintegration process was modelled as in [5]: for each A the cross 
section associated to one- and two-nucleon emission <Ti(e') and 02 (e') are ap¬ 
proximated by Gaussians in the range e( nin < e' < 30 MeV; at higher energies, 
30 MeV < e' < 150 MeV, the total cross section er 3 is approximated as a con¬ 
stant with branching rations taken from a table. 

Finally, r p io n , CMB ( 24 , T, z) is computed by 

T pion,CMB(Ar, 2 ) = (l + 2 ) 3 Ar p T o 1 nCMB (proton,(l + 0 )r ,2 = O), (15) 

where rV, n CMB for protons at z = 0 is interpolated from a table whose val¬ 
ues were obtained by numerically integrating r pion,EBL is computed as a 

function of T and 2 via 2D interpolation from a table obtained by numerically 
integrating (fl 2 l) : and TjA is computed by numerically integrating (TH 21 ) when 
needed. 

When a particle interacts, we sample the type of interaction, the probability 
of each type being pj = 7 V _1 / T tot ■ If the interaction is photodisintegration, one 
of the three channels a ±, 02 and 03 is similarly selected and, if the last channel 
is selected, the number of nucleons ejected is sampled from the branching ratios 
table. 
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Sampling number, type and energy of secondary particles 

Photodisintegration. When a nucleus is photodisintegrated, its energy 
is assumed to be split among the residual nucleus and the liberated nucleons 
proportionally to their mass, i.e. all the fragments inherit the Lorentz factor of 
the original nucleus. It is assumed that each nucleon has the same probability 
of being ejected, regardless of its type (i.e., if a nucleus with 26 protons and 30 
neutrons loses a nucleon, it is assumed to be a proton with probability 26/56 
and a neutron with probability 30/56); this simplifying assumption is only ap¬ 
proximately realistic (as in reality interaction channels yielding stable nuclei 
are more likely) and may result in a slight overestimate of the number of beta 
decays (and resulting neutrinos). 

As an exception, when a 9 Be nucleus with energy E is disintegrated, the 
fragments are two 4 He nuclei with energy 4E/9 each and a neutron with en¬ 
ergy E/0. Nuclei with mass numbers between 5 and 8 are never produced in 
our simulations. 


Pion photoproduction. When a pion is photoproduced, if the incoming 
particle is a nucleus, the nucleon that undergoes the interaction is chosen at 
random. We approximate all photo-hadronic processes as single-pion produc¬ 
tion; assuming isospin invariance, a neutral pion is produced (p + 7 —> p + n 0 
or n + 7 —> 7 T + 7 T 0 ) with probability 1/3 and a charged pion is produced 
(p + 7 —> n + n + orn + 7 —> p + n~) with probability 2/3. 

In order to sample the pion energy, first the photon energy e in the lab frame 
is sampled from its marginal distributiorQ 

p(0 de = + 4mFe)^de, ^ < e < +00 , (16) 

where m is the nucleon mass and eC n = m w + m 1 2 /2nv, then the squared centre- 
of-mass (CoM) energy s is sampled from its conditional distribution given e 


p(s|e) ds 


(s — m 2 )a(s) ds 

$(m 2 + 4mre) ’ 


(m + uItt) 2 < s < m 2 + AmTe , 


(17) 


from s the pion energy and momentum in the (CoM) frame are calculated as 


s-m 2 + ml * y/{s - (m + m J 2 ) (s - (m - m^) 2 ) 

K = —vi— 5 = -vi- (18) 

and the Lorentz factor of the transformation from the CoM frame to the lab 
frame as 7 = mV / y/s; then the pion energy is converted to the lab frame 
as El,. = 7 [E* +pV 0 S ^ 7 t)j where the distribution of 8„ is approximated as 
isotropic (cos0 w uniformly distributed between —1 and 1). In the lab frame, 
the momentum component orthogonal to the original travel direction is much 
smaller than that parallel to it, by a factor of order e/E ~ ICC 20 , so we neglect 
transverse components keeping the approach based on one dimensional propaga¬ 
tion. 

1 In practice, we use the fact that the marginal distribution of e corresponds to a distribution 

of 1(e) proportional to $(m 2 +4mre), and the conditional distribution of s given e corresponds 
to a uniform distribution of <&(s), so we actually sample I and $ and invert the functions to 
find the corresponding e and s. 
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The pion with energy E n , the nucleon with energy mT — E n , and (in the 
case of nuclei) a nucleus with mass number A — 1 and energy (A — l)mT are 
then added to the stack. 

3.2 Decay of unstable particles 

When an unstable particle is produced, it is assumed to decay instantaneously, 
as decay lengths are generally much shorter than all other relevant length scales. 
The energies of the decay products are sampled as described below and the decay 
products are added to the stack. 

Beta decay of neutrons and unstable nuclei Neutrons and nuclei not in 
the list of beta-decay stable isobars are assumed to immediately undergo beta 
decay. The Q-value of the reaction is read from a table taken from Ref. m or, 
for nuclei not on that table, estimated via the semi-empirical mass formula. 

The electron energy in the nucleus rest frame E* is sampled from a distri¬ 
bution oc (E* 2 — m 2 yi 2 El{Q — (E* — m e )) 2 (i.e., neglecting electromagnetic 
effects) and the neutrino energy is calculated as E* = Q — (E* — m e ); the recoil 
of the nucleus is neglected. The neutrino energy is converted to the lab frame by 
E u =TE*{ 1 — cos 0), where T is the Lorentz factor of the nucleus and 1 — cos 9 
is sampled from the uniform distribution between 0 and 2. 

The daughter nucleus (with the same energy and mass number A as the 
parent, with electric charge Z incremented in /3 - decay and decremented in 
/3 + decay) and the neutrino (D e in f3~ decay, v e in f3 + decay) are then added to 
the stack. 

Neutral pion decay A tt° with energy E„ decays into two photons with 
energy l? 7l distributed uniformly from 0 to E n and E l2 = E„ — E lx . 

Charged pion decay A ^ with energy E n decays into a muon with en¬ 
ergy E^ distributed uniformly from 0 to (1 — m^/m^.)^ and a neutrino with 
energy E v = E n - E^. 


Muon decay A muon with energy E ^ decays into two neutrinos and an elec¬ 
tron (ignored in SimProp); the energies E Vl ,E V2 of neutrinos are sampled as 
follows: 

• the energies of the neutrinos in the muon rest frame E* r , E* 2 are sampled 
independently uniformly from 0 to m^/2 — m 2 /2m IJ ,, and that of the elec¬ 
tron is E* =771^- E* ± - E* 2 ; 

• the corresponding momenta are computed as p* = E* t , p* 2 = E* 2 , and 

Pe = V E f - m e! 

• if these values violate any of the constraints E* > m e , p* Vl < p* 2 +Pe> 
Pt 2 < P* e + P* Vl , or p* e < p* Ul +pl 2 , they are discarded and a new E* ± , E* 2 
pair is sampled; 
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• the angle 8 i 2 between the two neutrinos is given by 


COS 012 = 


■ p * 2 


■ p * 2 


Zptipl 


(19) 


• the angle 8 1 between the first neutrino and the line of sight is isotropic, 
i.e. cos uniform from —1 to 1 ; 

• the angle </> between the second neutrino and the plane containing the line 
of sight and the first neutrino is uniform from 0 to 27r; 

• the angle 82 between the second neutrino and the line of sight is given by 
cos 82 = cos 0 i 2 cos 0 i — sin 0 i2 sin 0 i cos </>; 

• finally, neutrino energies are transformed to the lab frame via E Ul = 
7 {Et, + P*, cos 0i) and E V2 = 7 (E * 2 + p * 2 cos 0 2 ). 

3.3 Other particles: photons, electrons and neutrinos 

The propagation of photons and electrons produced is not yet implemented in 
SimProp: photons have their production energy and redshift recorded in the 
output file and electrons are neglected altogether; these studies will be matter 
of a forthcoming version of the SimProp code. 

Propagation of neutrinos is trivial: only the adiabatic energy loss due to 
the expansion of the Universe plays a role, being the opacity of the Universe 
to neutrinos irrelevant at red-shift 0 < 10 HU. Therefore, a neutrino produced 
with energy E at redshift z will reach Earth with energy E/{ 1 + 2 ). 

4 Structure of the output file 

The SimProp output is organized in ROOT JT2] files. In addition to the output 
tree of SimProp v2r0, SimProp v2r2 can optionally also write a tree named 
summary with the following branches: 

event A progressive number for the event, starting from 0. 

injEnergy The injection energy of the primary nucleus, in eV. 

injRedshift The redshift of the source. 

injDist The light-travel distance from the source to Earth, in Mpc. 

injA The mass number of the primary nucleus. 

injZ The atomic number of the primary nucleus. 

nNuc The total number of nuclei reaching Earth. 

nucEnergy Vector with the energies of the nuclei at Earth, in eV. 

nucA Vector with the mass numbers of the nuclei reaching Earth. 

nucZ Vector with the atomic numbers of the nuclei reaching Earth. 
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Figure 1: Fluxes of neutrinos at Earth expected in the Engel-Seckel-Stanev [13] 
scenario. The histograms were obtained via SimProp and the smooth curves 
refer to the data from Ref. m Left panel: v e (solid) and v e (dashed); right 
panel: v^ (solid) and (dashed). 


nPho The total number of photons produced by neutral pion decay. 

phoEProd Vector with the production energies of photons. 

phozProd Vector with the redshifts of the production points of photons. 

nNeu The total number of neutrinos and antineutrinos reaching Earth. 

neuEnergy Vector with the energies of neutrinos at Earth, in eV. 

neuFlav Vector with the flavours of neutrinos at production: +1 for u e , —1 for 
77 e , +2 for i/ , and —2 for (Neutrino oscillations are not implemented.) 


A Comparison with Engel—Seckel—Stanev fluxes 

In order to test our algorithm, particularly the new module of the code developed 
to compute the flux of secondary neutrinos, we have compared our results with 
those of a simulation performed by Engel, Seckel and Stanev m- The assump¬ 
tions of m are: (i) pure proton injection with spectrum oc E 2 exp(—E/10 21 ' 5 
eV) from 10 19 eV to 10 23 eV and total power density Pq = 4.5x 10 44 erg/Mpc 3 /yr 
at z = 0; (ii) no EBL; (iii) density of sources proportional to (1 + z) 3 for 
z < 1.9, constant for 1.9 < z < 2.7, and proportional to exp(— z/2.7) for 
2.7 < z < 8 . In figure [T] we plot the fluxes of secondary neutrinos and anti¬ 
neutrinos v e , v e , , results of [T3] are plotted as smooth lines while our 

results through histograms. 

The two simulations show results in quite good agreement. The only sub¬ 
stantial difference is that our electron antineutrino spectrum lacks the second 
peak at around 10 18 eV. This follows from our choice of neglecting the sub¬ 
dominant interaction channels in which several pions and/or heavier mesons 
are produced. This simplification is observationally irrelevant, because the phe¬ 
nomenon of neutrino oscillations (not shown here) will mix the flavours so that 
at Earth all neutrino flavours will show equal fluxes and the total antineutrino 
production at around 10 18 eV is dominated by 
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